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Abstract 

We present the results of full new calculation of radiocarbon 14 C production in the Earth atmosphere, using a numerical 
Monte-Carlo model. We provide, for the first time, a tabulated 14 C yield function for the energy of primary cosmic 
ray particles ranging from 0.1 to 1000 GeV/nucleon. We have calculated the global production rate of I4 C, which 
is 1.64 and 1.88 at/cm 2 /sec for the modern time and for the pre-industrial epoch, respectively. This is close to the 
values obtained from the carbon cycle reservoir inventory. We argue that earlier models overestimated the global 14 C 
production rate because of outdated spectra of cosmic ray heavier nuclei. The mean contribution of solar energetic 
particles to the global 14 C is calculated as about 0.25% for the modern epoch. Our model provides a new tool to 
calculate the 14 C production in the Earth's atmosphere, which can be applied, e.g., to reconstructions of solar activity 
in the past. 
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1. Introduction 



Radiocarbon 14 C is a long-living (half-life about 5730 years) radioactive nuclide produced mostly by cosmic rays 
in the Earth's atmosphere. Soon after production, it gets oxidized to 14 C02 and in the gaseous form takes part in 



the complex global carbon c ycle 



in many applications (e.g., 



Bolin et al 



Dorman, 



2004 



1979). Radiocarbon is important not only because it is used for dating 



Kromer 



20091) . but als o because it forms a primary method of paleo 



reconstructions of solar activity on the millen nial time scales (e.g. 



1989 



Bard et al 



1997 



Muscheler et al 



Stuiver and Ouav 



1980; 



Stuiver and Braziunas, 



2007b . An essential part of the solar activity reconstruction from radiocarbon 



data is computation of 14 C produc tion by cosmic rays in the Earth's atmosphere. First such computat i ons were 



perfo rmed in the 1960-1970s (e.g.. 



Lingenfelter 



1963 



Lingenfelter and Ramatv 



1970; 



Light et al 



1973 



O'Brien, 



1979) and were based on simplified numerical or semi-empirical methods. Later, full Monte-Carlo simulations of 
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O'Brien (1979) and 



Masarik and Beer 



the cosmic -ray induced atmospheri c cascade had be en p erforme d (Masarik and Beer, 1999; Masarik and Been, 12009). 
Most of earlier models, including 



1999) deal with a prescribed functional 



shape of the galactic cosmic ray spectrum, which makes it impossible to be applied to other types of cosmic ray 
spectra, e.g., solar energetic particles, supernova explosions, etc. A flexible approach includes calculation of the yield 
function (the number of cosmogenic nuclei produced in the atmosphere by the primary cosmic rays of the given type 
with the fixed energy and unit intensity ou tside the atmosphere), wh i ch can be convoluted with any given energy 



spect: 



r um of the primary cosmic ray s (e.g 



2008 



Kovaltsov and Usoskin, 



Webber and Higbie, 



2003 



Webber et al 



2007 



Usoskin and Kovaltsov 



20101) . This approach can be directly applied to, e.g., a problem of th e signatures of 



extreme solar energe t ic particle events i n the cosmogenic nucli de data, which is act i vely discussed (e.g., 



2006 



Hudson, 



2010; 



LaViolette, 



201 II) . Some earlier models (ILingen felter 



1963 



Castagnoli and Lai . 



Usoskin et al 



30) 



1980) provide 



the C yield function however it is limited in energy. Moreover, different models give results which differ by up to 
50% from each other, leading to large uncertainty in the global 14 C production rate. Therefore, the present status is 
that models providing the yield function are 30-50 years old and have large uncertainties. 

In addition, there is a systematic discrepancy between the results of theoretical models for the 14 C production 
and the global average 14 C production rate obtained from direct measurements of the specific I4 CC>2 activity in the 
atmosphere and from the carbon cycle reservoir inventory. While earlier production models predict that the global 



2s give systematically lower values ranging between 1 . 6 and 1 . 8 atoms/cm 2 /sec 



Damon and Sternberg, 



long (iLingenfeltei 



1989; 



O'Brien et al. 



1991 



1963) but is yet unresolved (iGoslai 



Goslar 



2001 



2001). 



Dorman, 



Lingenfelter 


1963 


Lai and Suess 


1968 



2004). This discrepancy is known since 



In this work we redo all the detailed Monte-Carlo computations of the cosmic-ray induced atmospheric cascade 
and the production of 14 C in the atmosphere to resolve the problems mentioned above. In Section|2]we describe the 
numerical model and calculation of the radiocarbon production. In Section [3] we compare the obtained results with 
earlier models. In Section [4] we apply the model to calculate the I4 C production by galactic cosmic rays and solar 
energetic particle events for the last solar cycle. Conclusions are presented in Section|5] 



ae 2. Calculation of the C production 

37 Energetic primary cosmic ray particles, when entering the atmosphere, collide with nuclei of the atmospheric 

38 gases initiating a complicated nucleonic cascade (also called shower). Here we are interested primarily in secondary 

39 neutrons whose distribution in the atmosphere varies with altitude, latitude, atmospheric state and solar activity. Neu- 

40 trons are produced in the atmosphere through multiple reactions including high-energy direct reactions, low-energy 

41 compound nucleus reactions and evaporation of neutrons from the final equilibrium state. Most of neutrons with en- 

42 ergy below 10 MeV are produced as an evaporation product of excited nuclei, while high-energy neutrons originate as 

43 knock-on neutrons in collisions or in charge exchange reactions of high-energy protons. While knock-on neutrons are 

2 



mainly emitted in the forward direction (viz. downwards), evaporated neutrons of lower energy are nearly isotropic. 
Radiocarbon 14 C is a by-product of the nucleonic cascade, with the main channel being through capture of secondary 
neutrons by nitrogen: N14(n,p)C14. Other channels (e.g., via spallation reactions) contribute negligibly, but are also 
considered here. 

We have performed a full Monte Carlo sim ulation of the nucleoni c component of the cosmic ray induced atmo- 



spheric cascade, using th e Planetocosmic co de (Desoreher 



particles through matter (!Gean t4 Collab oration et al 



et al 



20051) based on GEANT-4 toolkit for the passage of 



2003) (see details in Appendix). The secondary particles were 



tracked through the atmosphere until they undergo reactions with an air nucleus, exit the atmosphere or decay. In 
particular, secondary neutrons were traced down to epi-thermal energy. Simulations are computationally intensive. 
Simulations of single energies (ranging from 0.1 to 1000 GeV/nuc) were conducted, to determine the resulting flux of 
secondary neutrons. Since the calculations require very large computational time to keep the statistical significance 
of the results for low energies, we applied an analytical approach for atmospheric neutrons with energy below 10 
eV (see details in Appendix). Cross-sections have been adopted from the Experimental Nuclear Reaction Database 
(EXFOR/CSISRS) http://www.nndc.bnl.gov/exfor/exfor00.htm The number of simulated cascades induced by pri- 
mary CR particles was chosen as 10 5 - 10 6 to keep the statistical stability of the results at a reasonable computational 
time. Computations were carried out separately for primary protons and a-particles. Because of the similar rigid- 
ity /energy ratio, nuclei with Z > 2 were considered as effectively a-particles with the scaled number of nucleons (cf. 



Usoskin and Kovaltsov 



120081) . 

As the main result of these detailed computations we calculated the 14 C yield function. The yield functions for pri- 
mary protons and a-particles are tabulated in Table lA.fl and shown in Fig. lA.lK the energy range above 100 GeV/nuc 
is not shown). Note that the yields (per nucleon with the same energy) are identical for protons and cr-particle, viz. 
an a-particle is identical to four protons, at energies above 10 GeV/nuc. Details of the computations are given in 



Appendix Appendix A All further calculations are made using these yield functions. 

In order to compute the 14 C production q in the atmosphere at a certain place and conditions/time, one can use the 
following method: 

Xoo 
Y i (E)J i (E,t)dE, (1) 

where E is the particle's kinetic energy per nucleon, 7, is the spectrum of primary particles of type i on the top of the 
atmosphere, Ej C in GeV/nucleon is the kinetic energy per nucleon corresponding to to the local geomagnetic rigidity 
cutoff P c in GV. 

X , . , • (2) 



P c = -i ^E ic (E ic +2E r ), 

where E r = 0.938 GeV/nucleon is the proton's rest mass. Summation is over different types of the primary cosmic 
ray nuclei with charge Z, and mass A, numbers. T he local geomagnet ic rigidity cutoff is roughly defined via the 



geomagnetic latitude Aq of the location as following dElsasser et al 



1956) 



PJGV] = 1.9 • M • cos Ac, 



(3) 



where M is the dipole moment in units of [10 A m l of the Ear th's magnetic field. Although this approximation 



may slightly < 2% overestimate th e 14 C production (lO'Brien , 



(Dorman, 



2009; 



Clem et al 



2008), it is sufficient to study the global cosmic ray flux 



1997). The global production Q of radiocarbon is defined as the spatial global average 
of the local production q (both quantities give the number of 14 C nuclei produced per second per cm 2 of the Earth's 
surface). For the isotropic flux of primary particles in the interplanetary space (the level of anisotropy for galactic 
cosmic rays is usually smaller than 1%) the global production can be written as: 



poo 

«0= W 7. 



(E)J i (E,t)(l-f(E))dE, 



where the function 



RE) 



Jl - yfP(E)/ (1.9- M), ifP<1.9M 



(4) 



(5) 



0. 



if P > 1.9 M 



corresponds to sin(/lc) and accounts for the spatial average with the effect of the geomagnetic cutoff. 

Substituting any particular particle spectrum /, into Eq. |4] one can evaluate the 14 C production rate for different 
populations of cosmic rays, e.g., galactic cosmic rays (GCR), solar energetic particles (SEP), or more exotic sources 
like a nearby supernova explosion. 

First we consider the main source of 14 C, GCR modulated by the solar activity, using the standard approach. The 
energy spectrum of GCR particles of type i at 1 A U, ■/,-, is defined by the local interstellar spectrum (LIS), J\j[§ ,, and 



Usoskin et al 



the modulation potential (p as (see the formalism in 

J i (E,(f>) = JusAE + <I> i ) 



2005): 



(E)(E + 2E t ) 



(6) 



(E + ®i)(E + O, + 2E r ) 

where O,- = (eZ,/A,)</). The modulation potential cp is the variable related to solar activity, that parameterizes the shape 



' (e.g. 



Usoskin et al. 


2005 


Webber and Higbie 


20091 


Herbst et al. 


2010 


O'Brien 


2010) 



2010). Thus, the exact model of 



LIS must be specified tog ether with the values of <p. Here we use, as earlier, the proton LIS in the form (iBurger et al 



2000; 



Usoskin et al 



2005) 



Jus(E) 



1.9 x 10 4 -P(E) 



-2.78 



(7) 



1 + 0.4866 P(E)- 2 - 51 ' 

where P(E) = y/E(E + 2 E r ), J and E are expressed in units of particles/(m 2 sr s GeV/nucleon) and in GeV/nucleon, 
respectively. Here we consider two species of GCR separately: protons and heavier species, the latter including 
all particles with Z > 1 as a-particles with Z/A = 0.5 scaled by the number of nucleons. Heavier species should be 
treated separately as they are modulated in the heliosphere and Earth's magnetosphere differently, compared to protons 
because of the different Z/A ratio. Here we consider the nucleonic ratio of heavier particles (including a-particles) to 



2003 



Nakamura et al 



2010). 



protons in the interstellar medium as 0.3 (We bber and H igbie. 

The global 14 C production Q by GCR depends on two parameters, the solar magnetic activity quantified via the 
modulation potential <p and the Earth's geomagnetic field (its dipole moment M). The dependence is shown in the 



upper panel of Fig. [A3] One can see that both parameter s are equally important, and the knowledge of the geomagnetic 



field is very important (ISnowball and M uscheler 



20071) . In the lower panel, three cuts of the upper panel are shown 



to illustrate the effect of solar activity on Q, for the fixed geomagnetic field, corresponding to the modern conditions 
M - 7.8 ■ 10 22 A m 2 , as well as maximum (10 23 A m 2 ) and minimum (6 • 10 22 A m 2 ) dipole strength over the last 



ten millennia of the Holocene (Ko rte et al. 



201 1). The response of Q to changes of the geomagnetic field during the 



Holocene is within ±15%. However, the global I4 C would be nearly doubled during an inversion of the geomagnetic 
field ( viz. M — » 0). The modulation potential <p varies between about 300 and 1500 MV within a modern high solar 



2004 



cycle (Usoskin et al 



Usoskin et al 



2007 



201 1) , and can be as low as about 100 MV during the Maunder minimum (IMcCracken et al 



Steinhilber et al 



2008). Thus, changes of the solar modulation can also lead to a factor of 
2-3 variability on the global 14 C production rate. 

Next we investigated the sensitivity of Q to the energy of GCR. In Fig. IA.3I we show the relative cumulative 
production of 14 C, viz. the fraction of the total production caused by primary cosmic rays with energy below the given 
value E, as a function of E for different conditions. Often the median energy (the energy which halves the production) 



is used as a characteristic energy (e.g., 



Lockwood and Webber 



1996), which is the crossing of the curves in Fig. IA.3I 
with the horizontal dashed line. One can see that the median energy of 14 C production slightly changes with the level 
of solar activity, varying between 4 and 10 GeV/nuc corresponding to the Maunder minimum and the maximum of a 
strong solar cycle , respe ctively. The sensitivity of Q to the energy of GCR is close to that of a sea-level polar neutron 



monitor (cf. 



Beer 



2000). Slightly different shape of the neutron monitor cumulative response is due to the fact that it 



is ground-based while 14 C is produced in the entire atmosphere. 

As an example, we calculated the 14 C production predicted by the model for the last 60 years (see Fig . IA.4l i using 



122 the GCR modulation, reconstructed from the ground-based network of neutron monitors (lUsoskin et al 



2011), and 



IGRF (International Geomagnetic Reference Field - http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html) model of the 
Earth's magnetic field. The mean radiocarbon production for that period (1951-2010) is Q = 1.64 atom/cm 2 /sec, with 
the variability by a factor of two between 1.1 (in 1990 solar maximum) and 2.2 (in 2010 solar minimum) atom/cm 2 /sec. 
The me an 14 C production for the pre- industrial period (1750-190 0) calculated usin g the GCR modulation recon- 



struction by 



Alanko-Huotari et al 



Korte et al 



d201 II) is 1.88 atom/cm 2 /sec which 



(2007) and paleomagnetic data by 
is essentially lower than those reported in earlier works (1.9-2.5 atom/cm 2 /sec) and closer to the values obtained from 
the carbon cycle inventory (1 .6-1 .8 atom/c m 2 /sec) - see In troduction. This values can be further m 2% lower because 



of the used geomagnetic cut-off approach dO'Brien , 



2008). 



3. Comparison with earlier models 



In Fig. IA.ll we compare our present results with the yield functions calculated earlier (see the Figure caption for 
references). Our results are consistent with most of the earlier calculations (LR70 and DV91) within 10-20%. The 
CL80 yield function is not independently calculated but modified from LR70. While it is formally given for protons it 



effectively includes also ff-particles via scaling, thus being systematically higher than the other yield functions. Note 
that all the earlier computations of the yield function were limited in energy so that the upper considered energy of 
primary cosmic rays was from several to 50 GeV/nuc. On the other hand, contribution of higher energy cosmic rays 
is significant and may reach half of the total 14 C production (see Fig. [A3}. Here we present, for the first time, the 14 C 
yield function calculated up to TeV/nuc energy. Contribution from the higher energies is negligible because of the 
steep spectrum of GCR. 

Next we perform a more detailed comparison with the most recent 14 C production model by 



Masarik and Beer 



(120091 - MB09), who also used a GEANT-4 Monte-Carlo simulation tool. Since MB09 did not calculate the yield 
function, we use another way of comparison, via computing the global averaged 14 C production rate, as illustrated in 



Fig. IA.5I Our present result (black curve Q) in the Figure is systematical ly lower than that given by MB09 (big dots) 



145 by 25-30%. We suspect that the discrepancy arises fro m that 



for a prescribed GCR spectrum in the form given by (Garcia-Munoz et al 



Usoskin et al. 



1975 



2005 



Masarik and Beer (2009) calculated the 14 C pro duction 



Herbst et al 



Castagnoli and Lai, 



1980), which 



2010). Thus, in order to 



is different from the spectrum we use here (adopted from 
compare our results with those of MB09, we repeat our computations based on Eq.Q]to compute the global production 
Q* but using the same spectrum as MB09. The result for Q* is shown by the grey curves in Fig. IA.5I The overall 
agreement is within 5% but Q* value is systematically higher than that of MB09. The 5% difference can be related 
to the slightly different numerical scheme and also to the fact that MB09 treated an a-particle as four protons while 
we simulated them straightforwardly. In addition, the way of considering the geomagnetic shielding by MB09 is 
simplified (scaling) compared to our consideration (direct computations). We also compared the proton contributions 
(two dashed curves in Fig. IA.51 l to Q for the GCR spectrum discussed here (Eq. [6) and that used in MB09. The 
curves are nearly identical, suggesting that the difference in the used proton spectra is small and cannot be a cause 
for the observed systematic difference. We however notice a great difference between the a- (and heavier) particle 
spectra used here and in MB09. MB09 assumed of 12% for a— and 1% for heavier particle fr action in LIS (le ading 



to «0.64 nucleonic ratio between heavier species to protons in GCR) basing on the data from 



Simpson (1983). On 



the other hand, modern measurements (e.g., AMS, PAMELA) suggest that a-particles above 10 GeV/nuc contribute 
5-6% (in particle number) to LIS of GCR 
order of 0.25-0.3 outsid e the heliosphere (e.g 



Alcaraz et al 



2000a,b 



Nakamura et al 



... 



Adriani et al. 


2011 


Webber and Higbie, 


2003 



2010f) . viz. half of that assumed by MB09. Therefore, while we agree with MB09 in calculations 
of proton contribution into Q, they overestimate 14 C production by heavier species of GCR, using outdated spectra. 
This explains why the earlier results by MB99 and MB09 of 14 C production are systematically higher than our present 
result. 

Next we compare predictions of our model with other models' results for specific periods of time as shown in 
Fig. lA.6l (exact data sets used are mentioned in the Figure caption). One ca n see that our mo del predicts systematicall y 



O'Brienl (11979) and 



O'Brien et al 



(1991). 



lower production rates than most of other models, except of the model by j 
On the other hand, our yield function is generally consistent with others (Fig. IA.lt . indicating that the difference must 
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be related to the treatment of incoming GCR part icle spectra and /or geomagnetic shielding and not to the atmospheric 



cascade simulations. Models other than that by 



O'Brienl 0979) were based on theoretical calculations and included 



172 outdated overestimated abundance of a-particles, which explains the difference as discussed above. Therefore, we 

173 conclude that our model more correctly calculates the 14 C production as it agrees with the empirically-based models. 



74 4. 14 C production by solar energetic particles 

175 We also calculated prod uction of radiocarbon by solar energetic particles (SEP), because pres e ntly there is a wid e 



range of the results (e.g., 



Lingenfelter and Ramatv 



1970; 



Usoskinetal 



2006; 



Hudson, 



2010; 



LaViolette, 



2011). 



Here we compute the expected production of 14 C by the major known SEP events since 1951, using our calc ulated 



Tylka and Dietrich 



(2009). The 



yield function (Table I A. 1 b and SEP event-integrated spectra as reconstructed by 
corresponding production rate is shown by big open dots in Fig. lA.4l reduced to the monthly mean values. One can see 
that only a few SEP events can produce significant enhancements in 14 C production (« 70% in the monthly mean for 
the SEP event of 23-Feb-1956, 40% for 12-Nov-1960, 35% for two events in Oct-1989 and »20% for 29-Sep-1989). 
However, when applied to the annual time scale (the standard tree-ring time resolution), it gives only a few percent 
effect for years of maximum solar ac tivity and about 0.25% of the tota l contribution over the considered period. This is 



Lingenfelter and Ramatvl d 19701) (1.1% mean co ntribution of SEP into the global 



consistent with the earlier results by 
14 C production for 1954-1965, our model for the same period gives 0.8%) and by 



Usoskin et al 



(2006) (0.2% for 



1955-2005). Note that MB09, however, gives much smaller value of 0.02% for the SEP contribution to the global 
mean 14 C production, which i s probably caused by the neglect of the atmospheric cascade (and thus neutron capture 



188 channel) caused by SEPs (cf. 



Masarik and Reedv 



1995). 



5. Conclusions 



We have performed full new calculation, based on a detailed Monte-Carlo simulation of the atmospheric cascade 
by a GEANT-4 tool PLANETOCOSMIC, of the 14 C yield function. This is the first new calculation of the yield 
function since 1960-1970's, using modern techniques and methods, and the yield function is, for the first time 
ever, directly computed up to the energy of 1000 GeV/nuc (earlier models were limited to a few tens GeV/nuc 
and extrapolate d to higher ener gies). Our newly computed yield function gives the results which are in good 



agreement with 



O'Brien 



( 1979) and consistent with most of the earlier models, within 10-20%. 



We have calculated, using the new model and improved spectra of cosmic rays, the global production of 14 C, 
which appears to be significantly lower than earlier estimates and closer to the values obtained from the carbon 
cycle inventory. The calculated modern global production rate is 1.64 atom/cm 2 /sec, and the preindustrial rate 
(1750-1900 AD) is 1.88 atom/g/cm 2 , which is essentially lower than earlier estimates of 2-2.5 atom/cm 2 /sec. 
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We explain that the earlier models (including a recent model by iMasarik and Been (120091) ') overestimate the 
contribution of a-particle and heavier GCR species to the 14 C production, because of the use of outdated 
spectra. 

We have calculated, on the basis of the new m odel, contribution to the g lobal 14 C production by SEP events, 



using updated energy spectra reconstructions by 



Tylka and Dietrich 



(2009). The mean contribution of the SEPs 



for the last 50 years is estimated to be a;0.25% of the total production. 
• The present model provides an improved to ol to calculate the 14 C production in the Earth's atmosphere. Using 



the absolutely dated 14 C calibration curve (IReimer et al 



rays in the past (e.g., 



Solanki et al 



20091) . one can reconstruct the variability of cosmic 



2004) which, along with o ther long-term sol ar proxies has applications to 



paleoastrophysics, paleomagnetism and paleoclimatology (e.g.. 



Beer et al. 



2012). 
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Appendix A. Appendix: Details of numerical computations 



Numerical computations were done using the GEANT-based Monte-Carlo simulation tool Planetocosmic dDesorgher et al 



2005), which traces the atmospheric cascade induced by the primary cosmic ray partic les in full detail, inc luding the 



220 distribution of secondary particles. The Planetocosmic code has been recently verified (Usoskin et al 



2009) to agree 



1998), in the sense of en- 



within ss 10% with another commonly used Monte-Carlo package CORSIKA dHeck et al. 
ergy deposition in the atmosphere. The code simulates interactions and decays of various particles in the atmosphere 
in a wide range of energy. For the computations, we applied a realistic spherical atmospheric model NRMLSISE-00 



(Hedin, 



1991 



Picone et al 



2002). The QGSP_BIC_HP hadron interaction model has been applied with the standard 



electromagnetic interaction model. 

As an input for the simulations we used primary particles with fixed energy that impinge upon the top of the 
atmosphere at the random angle isotropically from the 2n solid angle. All computations were normalized per one 
such simulated particle. From the simulations we obtained the sum of secondary neutrons with energy within the AE 
energy bin centered at the energy E n , crossing a given horizontal level (atmospheric depth X g/cm 2 ), weighted with 
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230 1 1 / cos 6\ (where 6 is the zenith angle) to account for the geometrical factor, and divided by the energy bin width AE. 

231 This corresponds to the flux of secondary neutron with given energy F(E„,X) across a horizontal unit area, for the 

232 unit flux of primary cosmic rays on the top of the atmosphere. On the other hand, for quasi-stationary flux of neutrons 

233 this can be expressed as 

F(E„,X) = n„(E„,X)v„(E„), (Al) 

234 where n„ and v„ are the concentration (in [MeV cm 3 ] -1 ) and velocity of neutrons with energy E„ at the atmospheric 

235 depth level X. Let us denote the integral columnar flux as 

/(£„) = f " F(E n ,X)dX, (A2) 
Jo 

236 where X m = 1033 g/cm 2 is the total thickness of the atmosphere. Since our direct computations were performed down 

237 to energy of neutrons E\ =10 eV, we first computed the production of 14 C by these super-thermal neutrons, 

Gl = Z X(jC F(E "' X) nj(h) aj(En) dE "j dk ' (A3) 

238 where the outer integral is taken over the atmospheric height h, the concentration of target nuclei nfh) is defined as a 

239 product of the air density p and the content of the nuclei in a gram of air Kj, rtj(h) = p(h)Kfi o-j(E) is the cross-section 

240 of the corresponding reaction, and dX = p(h) dh, and summation is over target nuclei of different type (nitrogen 

241 kn = 3.225 • 10 22 atom/g; oxygen ko = 8.672 • 10 21 atom/g; argon kai = 1-94 • 10 20 atom/g, we also accounted for the 

242 isotopic distribution within these groups). Eqs. lA3l and lA2l can be transformed so that 

Gi = VjcJ I(E n )o-j(E n )dE n , (A4) 
J Jei 

243 All the cross-sections, used here, have been adopted from the Experimental Nuclear Reaction Database (EXFOR/CSISRS) 

244 http://www.nndc.bnl.gov/exfor/exforOO.htm 



245 We note that all the processes related to leakage of neutrons from the atmosphere (to the space or to soil) as well 

246 as their decay are accounted for in the direct simulation. 

247 Monte-Carlo simulations require extensive computational time in order to trace neutrons to thermal energy, thus 

248 compromising the statistical robustness of the results. On the other hand, the fate of 10 eV neutrons can be easily 

249 modelled theoretically, because of the simplicity of the processes involved, which allows us to save computational 

250 time and improve accuracy of the computations. The main process affecting epi-thermal neutrons in air is potential 

251 elastic scattering on N and O nuclei making neutrons to lose energy. After each elastic scattering, a neutron has a 

252 uniform distrib ution of ene rgy (in the laboratory frame) between its energy before the scattering E„ and a E„ (e.g. 

253 Chapter 7.2 in[ 



Fermi, 



2010). Here 



(A + l) 2 

254 where A is the mass number of the target nucleus. Then the probability for a neutron with the energy E„ (if E\ < 

255 E„ < E\la) before elastic scattering on a nuclei j to have energy E after the scattering so that E < E\ is {E\ - 
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256 ajE„)/(E„(l - aj)). Accordingly the "flux" (in the energy domain) of neutrons crossing the energy boundary E\ to 

257 (epi)thermal energies can be calculated as 

-Ei/aj £ _ a .£ 

. F(E n ,X)rij(X)a el j(E n ) 1 " 

!e, E„(l-aj) 



N=Y t f f ' ' F(E„,X) tij(X) cr el j(E n ) - ' ajE \ dE n dh (A6) 
^JhJE, E n (l-aj) 



258 or. using Eq. lA2l as 



ZrEi/aj £ _ a .£ 

kj I(E n )er el j(E„) 1 " dE„ (A7) 

Je, t ni.i -ap 



J Ei E„(l-aj) 

259 Reactions involving neutrons are: (1) N14(n,p)C14; (2) 017(n,a)C14; (3) N14(n,y)N15; (4) 016(n,y)017; (5) 

260 018(n,y)019 and (6) Ar40(n,y)Ar41. Note that only reactions (1) and (2) lead to production of 14 C while others 

261 simply provide a sink for neutrons. Cross-sections of neutron capture in all these reactions for energies below 10 eV 

262 can be expressed as 

crj = (A8) 

263 where Bj is a constant. Accordingly, the I4 C production by these neutrons can be calculated as 

G2=N *l-«N14 + *2-*017 (A9) 
I.j B J'*J 

264 The bulk of radiocarbon 14 C is produced via reaction (1) and about 0.001% in reaction (2). This is the main channel 

265 (95.8%) of the neutron sink. We have also considered leakage of neutrons from the upper atmospheric layers and decay 

266 of neutrons during their thermalization. These processes appear to be unimportant. In addition, we also computed 

267 possible contribution of secondary and primary protons to 14 C production via spallation reactions (e.g., 016(p,X)C14). 

268 These reactions are responsible for a negligible contribution to the total production. 

269 Then the final production of 14 C in the atmosphere by secondary neutrons corresponding to the primary cosmic 

270 ray particle with given energy is the sum of G\ and G2 and forms a point in the yield function Y/n. 
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Table A. 1 : Normalized yield functions Y p /n and Y a /n of the atmospheric columnar C production (in atoms sr) by a nucleon of primary cosmic 
protons and a-particles, respectively, with the energy given in GeV/nuc. For energy above 20 GeV/nuc, an a-particle is considered to be identical 
to four protons. 
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Figure A. 1: Yield function Y/n of 14 C production in the Earth's atmosphere by primary cosmic rays protons and a-particles (as denoted by 
"p" and "a" in th e legend, respectively) w ith given e nergy per nucleon. Different cu rves correspond to the pres ent work (T able lA.U and earlier 



models (CL80 ■ 
the legend. 



Castagnoli and La] , 



1980), (LR70 - 



Lingenfelter and Ramaty 



19701) and (DV91 - 



Dergach ev and Vekslei 



1991), as denoted in 
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Figure A.2: The global production of 14 C as function of the modulation potential <f> and the geomagnetic dipole moment M. The present value of 
M = 7.8 • 10 22 A m 2 is indicated by the thick arrow. The lower panel shows three cross-sections of the upper panel corresponding to the present 
value as well as to the maximum and minimum values of M over the past millennia, as indicated in the legend. Digital table for this plot is available 
at electronic supplement for this paper. 
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Figure A. 3: Relative cumulative production of C (fraction of the total production) as a function of the primary cosmic ray energy for different 
conditions: average solar activity (solid "aver." curve), solar maximum (tp = 1200 MV, dashed "S.max" curve), solar minimum (tp = 300 MV, 
dotted "S.min" curve), Maunder minimum (<p = 100 MV, circled "MM" curve). The thick grey curve corresponds to a polar sea-level neutron 
monitor. All curves are shown for the modern Earth magnetic field. 
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Figure A.4: Monthly averaged global production rates of 14 C since 1951 calculated using cosmic rays data from the world network of 
monitors and our calculated yield function. Open circles correspond to months with major solar energetic particle events. 
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Figure A. 5: Comparison of the global 14 C production rates, co mputed by different mod els as function of the modulation potential for the modern 

2009). Curves are computed using our calculated yield function 



Masarik and Beer 



geomagnetic field. Big dots correspond to the original results by 
(Table |A~"H and applying different cosmic rays spectra. Black curves (Q values) are calcu lated using the present results, while grey curves (Q* 
values) are calculated using our yield function but GCR spectra as used by I 



Masarik and Beer 



2009). Solid and dashed lines correspond to the total 



production and to production only by primary protons, respectively. 
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Figure A. 6: Global 14 C production as function of the modulation potential (f> as defined in 
th e present work's results. Symbols corresponds to ear lier works: 0'B79/91 (Ta b. 7 i n 



Lingenfelter 



1963 



Korff and Mendell. 



Lingenfelter and R amatv 



19801) ; MB09 (Tab. 3 in 



1970); Lal88 (Tabs. I and III in 



Masarik and Beeil 2009); MR95 (Tab. 1 in 



Usoskin e 



O'Brien, 



1979 



201 1). The thick grey curve presents 



1988); Lig73 (Tab. 6 in 



Masarik and Reedy, 1995). 



Light et al 



O'Brien et al., 1991); Lin63/LR70 (Tab. 1 



19731) . K80 (Tab. 1 in 



20 



